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SUMMARY 


The application of the vortex lattice method to problems not usually 
dealt with by this technique is considered. It is Bhown that if the dis- 
crete vortex lattice is considered as an approximation to surface-dis- 
tributed vorticity, then the concept of the generalized principal part 
of an integral yields a residual term to the vortex-induced velocity that 
renders the vortex lattice method valid for supersonic flow. Special 
schemes for simulating non-zero thickness lifting surfaces and fusiform 
bodies with vortex lattice elements are presented. Thickness effects of 
wing-like components are simulated by a double vortex lattice layer, and 
fusiform bodies are represented by a vortex grid arranged on a series of 
concentrical cylindrical surfaces. Numerical considerations peculiar 
to the application of these techniques are briefly discussed. 


INTRODUCTION 

The several versions or variations of the vortex lattice method that 
are presently available have proven to be very practical and versatile 
theoretical tools for the aerodynamic analysis and design of planar and non- 
planar configurations. The success of the method is due in great pert to 
the relative simplicity of the numerical techniques involved, and to the 
high accuracy, within the limitations of the basic theory, of the results 
obtained. But most of the work on vortex lattice methods appears to have 
concentrated on subsonic flow application. The applicability of the basic 
techniques of vortex lattice theory to supersonic flow has been largely 
ignored. It is one of the objectives of this paper to show how the vortex 
lattice method can be easily extended to deal with problems at supersonic 
Mach numbers with the same degree of success that it enjoys in subsonic flow. 

The other objective of this paper is to discuss a couple of schemes 
by which it is possible to simulate thickness and volume effects by using 
vortex lattice elements only. This represents an alternative, with somewhat 
reduced computational requirements, to the method of quadrilateral vortex 
rings (refs. 1 and 2). The simulation of thickness and volume effects makes 
possible the computation of the surface pressure distribution on wing-body 
configurations. The fact that this can be done without having to resort to 
additional types of singularities, such as sources, results in a simpler 
digital computer code. 


THE BASIC EQUATIONS 


Ward has shown, (ref. 3)» that the small-perturbation, linearized flow 
of an inviscid compressible fluid is governed by the three first order vector 
equations : 

V < v = F, V. w « Q, "w = T (l) 


on the assumption that the vorticity uT and the source intensity Q are known 
functions of the point whose position vector is R. The vector 'v is the 
perturbation velocity with orthogonal cartesian components u, v, and w, and 
'I' is a constant symmetrical tensor that for orthogonal cartesian coordinates 
with the x-axis aligned with the freestream direction has the form 
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where M w is the freestream Mach number. If 0 = 1-M ra , then the vector w 

has the components w’=@2ui+vj+wk. This vector was first introduced 

by Robinson (ref. 4), who called it the "reduced current velocity". If u 

denotes the total velocity vector, i.e. , u = (u» + u) T + v J + wl, and p 

the fluid density, then it can be shown that for irrotational and homentropic 

flow 


p u » p w u a + p a w + higher order terms 


(3) 


where the subscript ® indicates the value of the quantity at upstream infinity, 
e.g., uao 38 Ua, i. Therefore, to a linear approximation, the vector w is 
directly related to the perturbation mass flux as follows: 

w » ( p u - p„ u J /p. (4) 

The second equation of (l),. i.e., the continuity condition, shows that for 
source-free flows (Q * 0), w is a conserved quantity. 

Ward has integrated the three first order vector equations directly 
without having to resort to an auxiliary potential function. He obtained 
two different solutions for v (r), depending on whether 0^ is positive (sub- 
sonic flow), or negative (supersonic flow). These two solutions can be com- 
bined formally into a single expression if the following convention is used*. 

K * 2 for $ 2 > 0 

K = 1 for 0 2 < 0 X j 2 

Rp « Real part of |(x-X 2 _) 2 + 0 2 [(y-y 1 ) 2 + (z-z 1 ) 2 J| 


28 



I 


j . rr^sss ssmmmmmmmmm 


•f * Finite part of integral as defined by Hadamard (refs. 5 
and 6). 


The resulting solution for the perturbation velocity v at the point 
whose position vector is R^ * i + y^ j + z^ k, is given by 


v (E l ) ■ ‘ 2rtC 


-j~ n . w ( R) 


dS 
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dS 


ife «(«)V^dV^^^ixi(R)dV (5) 


This formula determines the value of v within the region V bounded by 
the surface S. The vector n is the unit outward (from the region V) normal 
to the surface S. Furthermore , it is understood that for supersonic flow 
only those parts of V and S lying within the domain of dependence (Mach 
forecone) of the point R^ are to be included in the integration. 


to 


mlm. 

For source-free (Q = 0), irrotational (uj=0) flow, equation (?) reduces 


(" 2 /' R-R 

< E l>"-2St ^-V.w(R)V|£ as + hjJ s ^ ldS(6> 


R-R, 

“ 

R- 


This is a relation between v inside S and the values of n.w and n >< v 
on S, but- these two quantities cannot be specified independently on S. 


To determine the source-fret:, irrotational flow about an arbitrary body 
B by means of equation (6), assume that the surface S coincides with the 
wetted surface of the body, with any trailing wake that it may have, and with 
a sphere of infinite radius enclosing the body and the whole flow field about 
it, namely, S * Sg + Ey + S . 


This surface S divides the 3pace into two regions, V e external to the 
body, and Vi internal to it. Applying equation (6) to both V e and Vi, since 
the integrals over S* converge to zero, the following expression is obtained* 
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vt.'Sj * 5 ^ k . i s (5) V;^ as - -^-)»x4iK5)|x -4p as (7) 

where ■ n^= -Kg is the unit normal to the body, or wake as the case may be, 
positive from the interior to the exterior of the body, & ^ " w e - wt, and 
4 v « v e - Here the subscripts designate the values of the quantities on 
the corresponding face of S. The first surface integral can be considered 
as representing the contribution of a source distribution of surface density 
N . A 7, while the second surface integral gives the contribution of a vorti- 
city distribution of surface density S x 4 ?. 

If the boundary condition of zero mass flux through the surface S_ + SL 
is applied to both external and internal flows 

N . p u e - ¥ . ( P(# i # + P w w e ) « 0 

N . p u t » N . (p a u # + « 0 ( 8 ) 

then the condition N . 4 w « 0 exists over + £L., and the flow field is 
uniquely determined by 

7(^5 “ - ^ dS (9) 

S B + S W 

where y (R) * N x 47 is the surface vorticity density. 


EXTENSION TO SUPERSONIC FLOW 


In order to extend the application of the vortex lattice method to 
supersonic flow, it is essential to consider the fundamental element of the 
method, the vortex filament, as a nv srical approximation scheme to the 
integral expression ( 9 ) instead of a real physical entity. The velocity field 
generated by a vortex filament can be obtained by a straightforward limiting 
process, the result being 


v(R x ) 


z£ 
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where f * 11 a Y • 6 

Y ^ oo 

6—o 

6 is a dimension normal to y> and di is the distance element along y. In the 
classical vortex lattice method, applicable only to subsonic flow, the vorti- 
city distribution over the body and the wake, i.e., over the surface 83 + %> 
is replaced by a suitable arrangement of vortex filaments whose velocity 
fields are everywhere determined by equation ( 10 ). This procedure is no 
longer appropriate for supersonic flow. For this latter case, it is necessary 
to go back to equation ( 9 ) and to derive an approximation to it. This is 
done in the following. 

If the surface S3 + fy, which defines the body and its wake, is considered 
as being composed of a large number of discrete flat area elements 7 over 
which the surface vorticity density Y can be assumed approximately constant, 
then equation (9) can be approximated by the following equation: 


N 



where N is the total number of discrete area elements T . When the point 
whose position vector is 81 is not part of Tj, the integral o/er this dis 
Crete area can be approximated by the mean value theorem as fellows: 



where Cj is a line in rj parallel to the average direction of Y in fj, 6 j 
is a distance normal to Cj, and di is the arc length element along Cj. Thip 
means that the velocity field induced by a discrete vorticity patch Tj can 
be approximated for points outside of t j by some mean discrete vortex line 
whov strength per unit length is yj V* But if the point Ri is part of the 
dl?cvete area T, the integral in equation (11) has an inherent singularity 
of the Cauchy type due to the fact that R ■ Ri at some point within r . In 
order to evaluate the Integral expression for this case, consider a point 
close to Ri but located Just above r by a distance e. As indicated in figure 
1 , the area of integration in 7 is divided into two regions, A r _ c and A e . 
Obviously, the integral over A T . e has no Cauchy-type singularity, Hadamard’s 
finite part concept being sufficient to perform the indicated integration. 
Thus, 
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The last integral in equation (13) represents the conventional discrete 
vortex line contribution whose evaluation presents no difficulty. In order 
to determine the integration denoted by I( ?) assume that, for simplicity, 
the coordinate system is centered at the point R]_, and that the x-y plane is 
determined by the discrete area r. Then, if y denotes the modulus of y. 


I( e) 




y sin A - x cos A 


!* 2 


- tf (y- 


+ e 


■r 


dx dy 


(14) 


where A is the angle between the y-axis and the direction of the vorticity 
in r , and * -g2 > 0 (supersonic flow). The components of the vector 
cross product y^(R-Ri) ■ y x R which are noc normal to the plane of r have 
been left out of equation (14) because, when the limit operation e-*o is 
carried out, they will vanish. The area Ac is bounded by a line parallel to 
the vorticity direction going through x -(1+B)e and by the intersection of the 
Mach forecone from the point (o> o, e) with the r -plane, consequently , if 
the integration with respect to x is performed first, 


l( e)=y cos A 
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-bJ/ 

f-* 


y + £ 


1x 2 -B 2 


(y 2 ♦ e 2 )f3/2 


dx | dy 


(15) 


ty -( 1+B) e 


where t = tan A , and Xi» X 2 are the values of y correspondi ng to the 
section of the line x=ty -(l+B)c with the hyperbola x _ _Wy2 + 2 . 

0 ■ € 2 (l-»2B)-2(l+B)cty -(B 2 -t 2 ) y 2 , then the finite part of the x- 

integration yields 


inter - 
Let 
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y cosA *X 


ty (ty-( 1+B) e) 

( y2 + ^2 ) >/^ 
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(16) 


Since e is a very small quantity, the variation of y in the interval ( Xj?) 
is going to be equally small, ana, therefore, the quantity within brackets in 
the last integrand of equation (l6) can be replaced by a mean value and taken 
outside of the integral sign. The same is not true of the term l/v0 since 
it will vary from oo for y = Xq, go through finite values in the integration 
interval, and then again increase to ao for y ® X2* With this in mind, and if 
y denotes a mean value of y, l( c) can be written as 

2 2 ~2 ^ 

X( g ) = v cos A B g-(l+B) cbv - (B -t ) y jT ( 1? ) 

B 2 y 2 + e 2 J Sf 


12 2 

But Xl , X2 are the roots of ty-e = -BVy + g , i.e., they are the roots of 
the polynomial denoted by <f>. Thus 

/T- \/c 2 (l + 2B)-2(l+B)Gty - (B 2 -t 2 ) y 2 =^B 2 -t 2 . \/ ( X^Ky-Xg) ( l8 ) 

Introducing this expression for J~<f into (17)» and taking the limit 
e-*-o, the following value for l( ") is obtained: 


l( o) = lim l( s) = - ^ C °p - — \j B 2 -t 2 -f— 

s — o B J 


_dy_ 


s/ ( X x - yXy-Xg) 


(19) 


The integral appearing in equation (19) can be easily evaluated by com- 
plex variable methods; its value is found to be 
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*2 


dy 


s /( \x-yXy-Xg) 


= TT 


( 20 ) 
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The contribution of the inherent singularity to the velocity field induced 
by vorticity patch T, vithin 7, denoted herein by w*, is therefore given by 


* 

v 


£ lim I( e) m m v cos A ,/J2 t 2 

2tt e -*o 2 v 


( 21 ) 


This contribution is perpendicular to the plane of T, and it has only 
physical meaning vhen b 2 > t 2 , i.e., vhen the vortex lines are swept in front 
of the Mich lines. It is expression (21), tak e n in conjunction with equation 
(12) , that makes the vortex lattice method applicable to supersonic flow. 


MODELING OF LIFTING SURFACES WITH THICKNESS 


The method of quadrilateral vortex rings placed on the actual body sur- 
face (ref. l) provides a way of computing the surface pressure distribution 
of arbitrary bodies using discrete vortex lines only. Numerical difficulties 
may occur when the above method is applied to the analysis of airfoils with 
sharp trailing edges due to the close proximity of two vortex surfaces of 
nearly parallel direction. An alternative approach, requiring somewhat less 
computer storage and easier to handle numerically, consists in using a double, 
or biplanar, sheet of swept horseshoe vortices to model a lifting surface 
with thickness, as shown schematically in figure 2. This constitutes an 
approximation to the true location of the singularities, similar in nature to 
the classical lifting surface theory approximation of a cambered sheet. 

All the swept horseshoe vortices, and their boundary condition control 
points, corresponding to a given surface, upper or lower, are located in a 
same plane. The upper and lower surface lattice planes are separated by a 
gap which represents the chordwise average of the airfoil thickness distri- 
bution. The results are not too sensitive to the magnitude of this gap; any 
value between one half to the full maximum chordwise thickness of the airfoil 
has been found to be adequate, the preferred value being two thirds of the 
maximum thickness. Furthermore, the gap can vary in the direction normal to 
the x-axis to allow for s panwise thickness taper. On the other hand, the 
chordwise distribution, or spacing, of the transverse elements of the horse- 
shoe vortices have a significant influence on the accuracy of the computed 
surface pressure distribution. For greater accuracy* for a given chordwise 
number of horseshoe vortices, the transverse legs have to be longitudinally 
spaced according to the 'cosine' distribution law 

l 1 - ("W)] (22 > 

where xj - x_ represents the distance from the leading edge to the midpoint 
of the swept leg of the Jth horseshoe vortex, c is the length of the local 
chord running through the midpoints of a given chordwise strip, and N is the 
number of horseshoe vortices per strip. The chordwise control point location 
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corresponding to this distribution of vortex elements is given by 



The control points are located along the center 1 no, or midpoint line, of the 
chordwise strip (fig. 3). Ian has shown (ref. 7 ) that the chordwise ’cosine' 
collocation of the lattice elements, defined by equations (22) and (23), 
greatly improve the accuracy of the computation of the effects due to lift. 
His results are directly extendable to the computation of surface pressure 
distributions of wings with thickness by the 'biplanar' lattice scheme pre- 
sented herein. 

The smalx perturbation boundary condition 

v . n' = -u .5 (24) 

00 

is applied at the control points. In equation (24), n * II + mj + nlc, and 
n' = mj + hk, where !, m, and n are the direction cosines of the normal to 
the actual airfoil surface. Equation (24) implies that l!u| « | mv + nw | . 
The use of the small perturbation boundary condition is consis .ent with the 
present 'biplanar' approach to the simulation of thick wings. 


MODELING OF FUSIFORM BODIES 


The modeling of fusiform bodies with horseshoe vortices requires a 
special concentrical vortex lattice if the simulation of the volume displace- 
ment effects, and the computation of the surface pressure distribution, are 
to be carried out. To define this lattice, it is necessary to consider first 
an auxiliary body, identical in cross-sectional shape and longitudinal area 
distribution to the actual body, with a straight barycentric line, i.e., 
without camber. The cross-sectional shape of this auxiliary body is then 
approximated by a polygon whose sides determine the transverse legs of the 
horseshoe vortices. The vertices of the polygon and the axis of the auxiliary 
body (which by definition is rectilinear (zero camber) and internal to all 
possible cross sections of the body) define a set of radial planes in which 
the bound trailing legs of the horseshoe vortices lie parallel to the axis 
(fig. 4). As the body cross section changes shape along its length, the 
corresponding polygon is allowed to change accordingly, but with the constraint 
that the polygonal vertices must always lie in the same set of radial planes. 
The axial spacing of the cross-sectional planes that determine the transverse 
vortex elements, or polygonal rings, follows the 'cosine' law of equation (22), 
The boundary condition control points are located on the auxiliary body sur- 
face, and in the bisector radial planes, with their longitudinal spacing given 
by equation (23). 

The boundary condition to be satisfied at these control points is the 
zero mass flux equation 


( 25 ) 


w . n ■ - u . n 

00 

2 

where all the components of the scalar product w.n = 0iu+mv+nw 
are to be retained. Thus, equation (25) is a higher order condition than 
equation (24). The use of this higher order boundary condition, within the 
framework of a linearized theory, is not mathematically consistent. There- 
fore, it can only be justified by its results rather than by a strict mathe- 
matical derivation. In the present treatment of fusiform bodies, it has been 
found that the use of higher order, or 'exact* boundary conditions is a re- 
quisite for the accurate determination of the surface pressure distribution. 

The fact that the vector w, instead of v, appears in the left hand member 
of equation (25) requires some elaboration. First, it should be pointed out 
that for small perturbations v . na v . n'. Furthermore, for incompressible 
flow (0 = l), the vector w is identical to the perturbation velocity v. Con- 
sequently, the boundary condition equation (24) is consistent with the con- 
tinuity equation, V. w = 0, to a first order f r compressible flow, and to any 
higher order for incompressible flow. But when a higher order boundary con- 
dition is applied in compressible flow to a linearized solution, it should be 
remembered that this solution satisfies the conservation of w, not of v, i.e., 
V. w = 0. Thus, the higher order boundary condition should involve the 
reduced current velocity, or perturbation mass flux, vector w, as in equation 
(25), rather than the perturbation velocity vector v. 

The body camber, which was eliminated in the definition of the auxiliary 
body, is taken into account in the computation of the direction cosines l,m, 
and n, which are implicit in equation ( 25 ). Therefore, the effect of camber 
is represented in the boundary condition but ignored in the spatial placement 
of the horseshoe elements. This scheme will give a fair approximation to 
cambered fusiform bodies provided that the amount of body camber is not too 
large. 


THE GENERALIZED VORTEX LATTICE METHOD 


Description of Method 

The three features discussed above, i.e., the inclusion of the vorticity- 
induced residual term w* for supersonic flow, the 'biplanar* scheme for rep- 
resenting thickness, and the use of a vortex grid of concentrical polygonal 
cylinders for the simulation of fusiform bodies, have been implemented in a 
computational procedure herein known as the Generalized Vortex Lattice (GVL) 
method. The GVL method has been codified in a Fortran IV computer program 
(VORIAX), which has been widely utilized throughout the Lockheed-California 
Company as an efficient aerodynamic design tool for advanced aircraft confi- 
gurations in subsonic and supersonic flows. 
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The basic element of the method is the swept horseshoe vortex with 
'bound' and 'free' legs. In the present version of the method, the free legs 
may trail to downstream infinity in any arbitrary, but predetermined, direc- 
tion. The lattice formed by the bound legs of the horseshoe vortices is 
laid out on the proper cylindrical surfaces, the trailing legs being parallel 
to the x-axis. Figure 5 illustrates schematically the representation of a 
simple wing-body configuration within the context of the present method. The 
streamwise arrangement of the lattice follows the 'cosine' distribution law 
(eq. (22)), whereas the spanwise, or cross-flow, spacing of the trailing legs 
can be arbitrarily specified. To each horseshoe vortex there corresponds an 
associated control point, placed midway between the bound trailing legs of 
the horseshoe and longitudinally spaced according to equation (23). 

The velocity field induced by the elementary horseshoe vortex is derived 
from equation (12), and it includes the contribution given by equation (21) 
when the velocity induced by a horseshoe at its own control point is evaluated 
at supersonic Mach numbers. This velocity field is used to generate the co- 
efficients of a system of linear equations relating the unknown vortex 
strengths to the appropriate boundary condition at the control points. This 
linear system is solved by either a Gauss-Seidel iterative procedure (ref. 8), 
or by a vector orthogonalization technique (ref. 9). 

The pressure coefficients are computed in terms of the perturbation 
velocity components. Force and moment coefficients are determined through a 
numerical integration process. Due account is taken of the leading edge 
suction through the application of Lan's procedure (ref. 7), which the GVL 
method directly extends to supersonic flow. 


Numerical Considerations 

At supersonic Mach humbers , the velocity induced by a discrete horseshoe 
vortex becomes very large in the very close proximity of the envelope of Mach 
cones generated by the transverse leg of the horseshoe. At the characteristic 
envelope surface itself, the induced velocity correctly vanishes, due to the 
finite part concept. This singular behavior of the velocity field occurs only 
for field points off the plane of the horseshoe. For the planar case, the 
velocity field is well behaved in the vicinity of the characteristic surface. 

A simple procedure to treat this numerical singularity consists of defining 
the characteristic surfaces by the equation 

(x-x ^ 2 = c B 2 | (y-yj ) 2 + (z-zi) 2 | ( 26 ) 

where C is a numerical constant whose value is greater than, but close to, 1. 
It has been found that this procedure yields satisfactory results, and that 
these results are quite insensitive to reasonable variations of the parameter 

C. 
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Another numerical problem, peculiar to the supersonic horseshoe vortex, 
exists in the planar case (field point in the plane of the horseshoe) when 
the field point is close to a transverse vortex leg swept exactly parallel 
to the Mach lines (sonic vortex), while the vortex lines immediately in front 
of and behind this sonic vortex are subsonic and supersonic, respectively. 
This problem can be handled by replacing the boundary condition equation for 
such sonic vortex with the averaging equation 

+ 2 Yj* " Yj* +1 ■ 0 (27) 


where y T * is the circulation strength of the critical horseshoe vortex, and 
Y I *_ 1 an ^ Yj* +1 are the respective circulation values for the fore-and-aft 

adjacent subsonic and supersonic vortices. 


The axialwash induced velocity component (u) is needed for the computa- 
tion of the surface pressure distribution, and for the formulation of the 
boundary condition for fusiform bodies. When the field point is not too 
close to the generating vorticity element, the axialwash is adequately des- 
cribed by the conventional discrete horseshoe vortex representation. But if 
this point is in the close vicinity of the generating element, as may occur 
in the biplanar and in the concentrical cylindrical lattices of the present 
method, the error in the computation of the axialwash due to the discretiza- 
tion of the vorticity becomes unacceptable. This problem is solved by resor- 
ting to a vortex-splitting technique, similar to the one presented in refer- 
ence 10. Briefly, this technique consists of computing the axialwash induced 
by the transverse leg of a horseshoe as the summation of several transverse 
legs longitudinally redistributed, according to an interdigitation scheme, 
over the region that contains the vorticity represented by the single discrete 
vortex. This is done only if the point at which the axialwash value is re- 
quired lies within a given near field region surrounding the original dis- 
crete vortex. 
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COMPARISON WITH OTHER THEORIES AND EXPERIMENTAL RESULTS 
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Conical flow theory provides a body of 'exact* results, within the con- 
text of linearized supersonic flow, for some simple three-dimensional confi- 
gurations. These exact results can be used as bench mark cases to evaluate 
the accuracy of numerical techniques. This has been done rather extensively 
for the GVL method, and very good agreement between it and conical flow theory 
has been observed in the computed aerodynamic load distribution and all force 
and moment coefficients. Only some typical comparisons are presented in this 
paper, figures 6 through 9* 

Finally, the capability of computing surface pressure distributions by 
the method of this paper is illustrated in figures 10 and 11. 
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CONCLUDING REMARKS 






It has been shown that vortex lattice theory can be extended to super- 
sonic flow if due account is taken of the principal, part of the surface vorti- 
city integral. Furthermore , special vortex lattice layouts, which allow the 
simulation of thickness and volume with horseshoe vortices, have been presen- 
ted. All this greatly enhances the value of vortex lattice theory as a com- 
putationally efficient design and analysis tool, as exemplified by its exten- 
sive vise at the Lockheed-California Company, discussion of which has been 
precluded by space limitations. 
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Figure 2.- Modeling of thick wing with horseshoe vortices. 
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Figure 3.- Vortex lattice collocation. 
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AXIAL DISTRIBUTION OF B.C. CONTROL POINTS: 

Figure 4.- Modeling of fusiform body with horseshoe vortices 
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Figure 6.- Theoretical comparison of arrow wing lift slope 
and aerodynamic center location. 
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Figure 7.- Theoretical comparison of arrow wing drag-due-to-lif t factor. 
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Figure 8.- Theoretical comparison of chordwise 
loading for delta wing. 
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Figure 9.- Theoretical comparison of chordwise loading 
for sweptback rectangular wing. 
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